Analysis for ‘Functional traits predict species responses to environmental variation in a California grassland annual plant community’
Gaurav Kandlikar
12/21/2021
requirements <- c('here', 'bbmle', 'broom', 'broom.mixed',
'DHARMa', 'effects', 'factoextra', 'FactoMineR', 'TMB',
'glmmTMB', 'latex2exp', 'patchwork',
'tidyverse', 'visreg', 'devtools', 'performance', 'kableExtra')
lapply(requirements,function(x){if(!require(x,character.only = T)){install.packages(x)}})
# Install data package for cover data:
devtools::install_github("akleinhesselink/sedgwickcover")
knitr::opts_chunk$set(echo = TRUE, warning = FALSE)
knitr::opts_knit$set(root.dir = here::here())R Code to recreate all analyses in the manuscript
Setup
library(bbmle)
library(broom)
library(broom.mixed)
library(DHARMa)
library(effects)
library(factoextra)
library(FactoMineR)
library(glmmTMB)
library(latex2exp)
library(patchwork)
library(tidyverse)
library(visreg)
library(paran)
library(sedgwickcover)
theme_plots <- function() {
theme_minimal() +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
plot.tag = element_text(face = "bold"),
axis.title = element_text(size = 12),
plot.caption = element_text(size = 12)
)
}
theme_set(theme_plots())Community-weighted trait means
Import trait and cover data
strts_all <- read_csv("data/all_trait_data.csv") %>%
dplyr::select(-`...1`)## Rows: 65 Columns: 15## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): USDA_symbol, dataset
## dbl (13): ...1, leaf_size, SLA, LDMC, seed_mass, max_height, SRL, relative_s...##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.strts_all %>%
group_by(USDA_symbol) %>%
summarize(n = n()) %>%
arrange(-n)# Remove some duplicate rows from the traits df
# (These species were measured in multiple years;
# retaining only those measurements made in 2017, i.e.
# the year of the community cover sampling).
strts <- strts_all %>%
group_by(USDA_symbol) %>%
filter(!(USDA_symbol == "ANAR" & dataset == "TAPIOCA")) %>%
filter(!(USDA_symbol == "BRMA3" & dataset == "2016")) %>%
filter(!(USDA_symbol == "CHGL" & dataset == "2016")) %>%
filter(!(USDA_symbol == "FEMI2" & dataset == "2016")) %>%
filter(!(USDA_symbol == "HOMU" & dataset == "2016")) %>%
filter(!(USDA_symbol == "MICA" & dataset == "2016")) %>%
filter(!(USDA_symbol == "MIDO" & dataset == "TAPIOCA")) %>%
filter(!(USDA_symbol == "SIGA" & dataset == "TAPIOCA")) %>%
return
# Import the cover data
p17 <- read_csv("data/sedgwick_cover_2017.csv")## Rows: 4400 Columns: 6## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): USDA_symbol
## dbl (5): year, site, plot, area, cover##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Calculate the relative cover in each plot at each site
p17 <- p17 %>%
group_by(site, plot) %>%
summarize(total_cover = sum(cover)) %>%
left_join(p17, .) %>%
mutate(relative_cover = cover/total_cover) %>%
as_tibble()## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.## Joining, by = c("site", "plot")# Merge the cover and trait data frames
p17_trt <- left_join(p17, strts)## Joining, by = "USDA_symbol"# species with missing data for at least one of the traits
missing_trts <-
p17_trt %>%
group_by(USDA_symbol, site, plot) %>%
filter(is.na(SLA) | is.na(SRL) | is.na(max_height)| is.na(seed_mass)) %>%
arrange(-relative_cover)
missing_trts_nz <- missing_trts %>% filter(relative_cover > 0)
mean(missing_trts_nz$relative_cover)## [1] 0.04041392# std error function
sem <- function(x) (sd(x)/sqrt(n()))
# Summarize the trait data to begin calculating CWM mean
# (and sd of CWM trait) at each site
p17_trt_sum <-
p17_trt %>%
mutate(lsla = (log(SLA)),
lsrl = (log(SRL)),
lmaxht = (log(max_height)),
lseedmass = (log(seed_mass))) %>%
dplyr::select(site:relative_cover, lsla, lsrl, lmaxht, lseedmass) %>%
mutate(lsla_contrib = lsla*relative_cover,
lsrl_contrib = lsrl*relative_cover,
lmaxht_contrib = lmaxht*relative_cover,
lseedmass_contrib = lseedmass*relative_cover) %>%
group_by(site, plot) %>%
summarize_at(c("lsla_contrib", "lsrl_contrib", "lmaxht_contrib", "lseedmass_contrib"),
sum, na.rm = T) %>%
rename(lsla_cwm = lsla_contrib,
lsrl_cwm = lsrl_contrib,
lmaxht_cwm = lmaxht_contrib,
lseedmass_cwm = lseedmass_contrib)
# Summarize across the sub-plots at each site, to get the
# mean CWM, and SD+SEM at each site
p17_sitesum <-
p17_trt_sum %>%
group_by(site) %>%
summarize_at(c("lsla_cwm", "lsrl_cwm", "lmaxht_cwm", "lseedmass_cwm"),
lst(mean, sd, sem))
env <- read_csv("data/all_environmental_data.csv")## Rows: 16 Columns: 22## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): type, microsite
## dbl (20): plot, lat, lon, ele, Tmax, Tmin, organic_matter_ENR, pH, CEC_meq_1...##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.env$ca_mg <- env$Ca_ppm/env$Mg_ppm
# merge in the environmental data with both plot-level and site-level CWM dfs
p17_sitesum_env <- env %>%
rename(site = plot) %>%
left_join(p17_sitesum, .)## Joining, by = "site"Now, compute CWM-environment relations:
# Now, test for CWM trait shifts along the focal environmental axes ----
# Import functions that run LMs and make plots
source("functions/make_trait_env_lms.R")
source("functions/make_trait_env_plot.R")
lsla_sand_lm <- make_trait_env_lms(trait = "lsla_cwm_mean", env = "sand")
lsla_sand_plot <- make_trait_env_plot(trait = "lsla_cwm_mean", env = "sand",
lsla_sand_lm, xl = "Sand (%)", yl = "CWM log(SLA (cm^2/g))")
lsla_camg_lm <- make_trait_env_lms(trait = "lsla_cwm_mean", env = "ca_mg")
lsla_camg_plot <- make_trait_env_plot(trait = "lsla_cwm_mean", env = "ca_mg",
lsla_camg_lm, xl = "Ca:Mg", yl = "CWM log(SLA (cm^2/g))")
lsla_depth_lm <- make_trait_env_lms(trait = "lsla_cwm_mean", env = "depth")
lsla_depth_plot <- make_trait_env_plot(trait = "lsla_cwm_mean", env = "depth",
lsla_depth_lm, xl = "Depth (cm)", yl = "CWM log(SLA (cm^2/g))")
lsrl_sand_lm <- make_trait_env_lms(trait = "lsrl_cwm_mean", env = "sand")
lsrl_sand_plot <- make_trait_env_plot(trait = "lsrl_cwm_mean", env = "sand",
lsrl_sand_lm, xl = "Sand (%)", yl = "CWM log(SRL (m/g))")
lsrl_camg_lm <- make_trait_env_lms(trait = "lsrl_cwm_mean", env = "ca_mg")
lsrl_camg_plot <- make_trait_env_plot(trait = "lsrl_cwm_mean", env = "ca_mg",
lsrl_camg_lm, xl = "Ca:Mg", yl = "CWM log(SRL (m/g))")
lsrl_depth_lm <- make_trait_env_lms(trait = "lsrl_cwm_mean", env = "depth")
lsrl_depth_plot <- make_trait_env_plot(trait = "lsrl_cwm_mean", env = "depth",
lsrl_depth_lm, xl = "Depth (cm)", yl = "CWM log(SRL (m/g))")
lmaxht_sand_lm <- make_trait_env_lms(trait = "lmaxht_cwm_mean", env = "sand")
lmaxht_sand_plot <- make_trait_env_plot(trait = "lmaxht_cwm_mean", env = "sand",
lmaxht_sand_lm, xl = "Sand (%)", yl = "CWM log(Max. Height (cm))")
lmaxht_camg_lm <- make_trait_env_lms(trait = "lmaxht_cwm_mean", env = "ca_mg")
lmaxht_camg_plot <- make_trait_env_plot(trait = "lmaxht_cwm_mean", env = "ca_mg",
lmaxht_camg_lm, xl = "Ca:Mg", yl = "CWM log(Max. Height (cm))")
lmaxht_depth_lm <- make_trait_env_lms(trait = "lmaxht_cwm_mean", env = "depth")
lmaxht_depth_plot <- make_trait_env_plot(trait = "lmaxht_cwm_mean", env = "depth",
lmaxht_depth_lm, xl = "Depth (cm)", yl = "CWM log(Max. Height (cm))")
lseedmass_sand_lm <- make_trait_env_lms(trait = "lseedmass_cwm_mean", env = "sand")
lseedmass_sand_plot <- make_trait_env_plot(trait = "lseedmass_cwm_mean", env = "sand",
lseedmass_sand_lm, xl = "Sand (%)", yl = "CWM log(seed mass g)")
lseedmass_camg_lm <- make_trait_env_lms(trait = "lseedmass_cwm_mean", env = "ca_mg")
lseedmass_camg_plot <- make_trait_env_plot(trait = "lseedmass_cwm_mean", env = "ca_mg",
lseedmass_camg_lm, xl = "Ca:Mg", yl = "CWM log(seed mass g)")
lseedmass_depth_lm <- make_trait_env_lms(trait = "lseedmass_cwm_mean", env = "depth")
lseedmass_depth_plot <- make_trait_env_plot(trait = "lseedmass_cwm_mean", env = "depth",
lseedmass_depth_lm, xl = "Depth (cm)", yl = "CWM log(seed mass g)")
cwmtrts_plot <- wrap_plots(list(lsla_sand_plot, lsla_camg_plot, lsla_depth_plot,
lsrl_sand_plot, lsrl_camg_plot, lsrl_depth_plot,
lmaxht_sand_plot,lmaxht_camg_plot, lmaxht_depth_plot,
lseedmass_sand_plot, lseedmass_camg_plot, lseedmass_depth_plot),
nrow = 4)
cwmtrts_plot
cwmtrts_plot_seed <- wrap_plots(list(lseedmass_sand_plot, lseedmass_camg_plot, lseedmass_depth_plot),
nrow = 3)There are five significant CWM trait-environment relationships: SLA-Ca:Mg; SLA:Depth; SRL:Sand; Max. Height:Depth; Seed Mass:Depth
Demography experiment
Import performance data, merge with trait data, and merge with environmental variables data.
## Now, we can test for trait-environment interactions in the demography data ----
# Clean/manage performance data ------
performance <- read_csv("data/seed_production_processed_nocomp.csv")## Rows: 1399 Columns: 5## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): plot_type, sp_code
## dbl (3): plot_num, replicate, num_seeds_produced##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.performance <- performance %>%
# NA in this means the seed production was zero
mutate(num_seeds_produced = ifelse(is.na(num_seeds_produced), 0,
num_seeds_produced))
germination_data <- read_csv("data/germination_data.csv") %>%
rename(replicate = rep)## Rows: 1360 Columns: 6## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): sp_code
## dbl (5): plot_num, rep, num_germinants, expected_germinants_l, num_germinant...##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.performance <- left_join(performance, germination_data) ## Joining, by = c("plot_num", "replicate", "sp_code")performance <- performance %>%
# Make a new column for Hummock or matrix
mutate(plot_env = ifelse(plot_num %in% c(742:745,747,748,750,751,752,755),
"hummock", "matrix")) %>%
return
# Make a new df to help match species names between trait and performance dfs
sp_code_usda_symbol <- data.frame(
sp_code = unique(performance$sp_code),
USDA_symbol = c("URLI5", "AMME", "BRMA3", "CEME2",
"CLBO", "CLPU2", "CHGL", "EUSP",
"HECO7", "HOMU", "LACA7", "LOWR2",
"MEPO3", "MICA", "PLER3", "SACO6", "FEMI2")
)
performance <- left_join(performance, sp_code_usda_symbol)## Joining, by = "sp_code"# Add in the environmental data
perf_env <- left_join(performance,
env %>%
# keep only the focal environmental variables, for now
dplyr::select(plot_num = plot, unscaled_ca_mg = ca_mg,
unscaled_sand = sand, unscaled_depth = depth))## Joining, by = "plot_num"perf_env# Add in the trait data
# Include only the trait data collected in 2016 or for TAPIOCA
trait_df16 <-
strts_all %>%
filter(dataset %in% c("TAPIOCA", "2016")) %>%
dplyr::select(USDA_symbol, SLA, SRL, max_height, seed_mass)
# We only have seed mass from other years for BRMA and HOMU, so import those in.
trait_df16[which(trait_df16$USDA_symbol == "HOMU"), "seed_mass"] <-
strts_all[which(strts_all$USDA_symbol == "HOMU"), "seed_mass"] %>% unlist %>% .[1]
trait_df16[which(trait_df16$USDA_symbol == "BRMA3"), "seed_mass"] <-
strts_all[which(strts_all$USDA_symbol == "BRMA3"), "seed_mass"] %>% unlist %>% .[1]
perf_env_trt <- left_join(perf_env, trait_df16, by = "USDA_symbol")
perf_env_trt$unscaled_log_sla <- log(perf_env_trt$SLA)
perf_env_trt$unscaled_log_srl <- log(perf_env_trt$SRL)
perf_env_trt$unscaled_log_maxht <- log(perf_env_trt$max_height)
perf_env_trt$unscaled_log_seedmass <- log(perf_env_trt$seed_mass)
# perf_env_trt <- left_join(perf_env, traits)
perf_env_trt$plot_num <- as.factor(perf_env_trt$plot_num)
perf_env_trt$num_seeds_produced <- as.integer(perf_env_trt$num_seeds_produced)Now, run the trait-environment GLMMs
# Trait-Env model based on CWM trait ----
# SLA-Ca:Mg; SLA:Depth; SRL:Sand; Max. Height:Depth
# scale the relevant columns
perf_env_trt <- perf_env_trt %>% mutate(
ca_mg = scale(unscaled_ca_mg),
depth = scale(unscaled_depth),
sand = scale(unscaled_sand),
log_sla = scale(unscaled_log_sla),
log_srl = scale(unscaled_log_srl),
log_maxht = scale(unscaled_log_maxht),
log_seedmass = scale(unscaled_log_seedmass),
germ_ratio = num_germinants_corrected/expected_germinants_l,
germ_prod = floor(num_seeds_produced*germ_ratio))
trait_model_full <- glmmTMB(num_seeds_produced ~
ca_mg* log_sla +
depth*log_sla +
sand*log_sla +
ca_mg*log_srl +
depth*log_srl +
sand*log_srl +
ca_mg*log_maxht +
sand*log_maxht +
depth*log_maxht +
ca_mg*log_seedmass +
sand*log_seedmass +
depth*log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
family = nbinom2(link = "log"))
summary(trait_model_full)## Family: nbinom2 ( log )
## Formula:
## num_seeds_produced ~ ca_mg * log_sla + depth * log_sla + sand *
## log_sla + ca_mg * log_srl + depth * log_srl + sand * log_srl +
## ca_mg * log_maxht + sand * log_maxht + depth * log_maxht +
## ca_mg * log_seedmass + sand * log_seedmass + depth * log_seedmass +
## (1 | sp_code) + (1 | plot_num)
## Zero inflation: ~1
## Data: perf_env_trt
##
## AIC BIC logLik deviance df.resid
## 13246.9 13372.7 -6599.4 13198.9 1375
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## sp_code (Intercept) 0.4506 0.6713
## plot_num (Intercept) 0.2353 0.4851
## Number of obs: 1399, groups: sp_code, 17; plot_num, 16
##
## Dispersion parameter for nbinom2 family (): 1.54
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.27057 0.20738 25.415 < 2e-16 ***
## ca_mg 0.36426 0.15086 2.415 0.015756 *
## log_sla -0.18257 0.17028 -1.072 0.283640
## depth 0.32482 0.13719 2.368 0.017902 *
## sand -0.02416 0.13867 -0.174 0.861676
## log_srl 0.06550 0.18259 0.359 0.719806
## log_maxht 0.61972 0.18198 3.405 0.000661 ***
## log_seedmass 0.22919 0.17385 1.318 0.187384
## ca_mg:log_sla 0.13332 0.03418 3.901 9.6e-05 ***
## log_sla:depth -0.02694 0.02942 -0.915 0.359941
## log_sla:sand 0.04451 0.03097 1.437 0.150592
## ca_mg:log_srl -0.03663 0.04256 -0.861 0.389354
## depth:log_srl -0.01212 0.03705 -0.327 0.743616
## sand:log_srl -0.11098 0.03474 -3.195 0.001399 **
## ca_mg:log_maxht -0.04342 0.04026 -1.079 0.280739
## sand:log_maxht 0.05861 0.03865 1.516 0.129429
## depth:log_maxht 0.05954 0.04137 1.439 0.150109
## ca_mg:log_seedmass 0.07007 0.04136 1.694 0.090263 .
## sand:log_seedmass 0.01030 0.03517 0.293 0.769630
## depth:log_seedmass 0.08618 0.03603 2.392 0.016750 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.53232 0.05571 -9.556 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1trait_model_null <- glmmTMB(num_seeds_produced ~
1 +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
family = nbinom2(link = "log"))
trait_model_additive <- glmmTMB(num_seeds_produced ~
ca_mg + sand + depth +
log_sla + log_srl + log_maxht + log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
family = nbinom2(link = "log"))Check model performance:
# Verify whether interaction terms are supported
# relative to a model with only additive terms
# and a null model
bbmle::AICctab(trait_model_null, trait_model_full, trait_model_additive)## dAICc df
## trait_model_full 0.0 24
## trait_model_additive 21.7 12
## trait_model_null 30.9 5# Check VIF in the conditional component of the full model
performance::check_collinearity(trait_model_full, component = "conditional")# Check other model diagnostics using DHARMa
resids <- DHARMa::simulateResiduals(trait_model_full)
plot(resids)
# Looks some issue with uniformity - though the QQ looks fairly
# linear, and Florian Hartig's notes (e.g. https://github.com/florianhartig/DHARMa/issues/181)
# suggest that significance of Kolmogorov-Smirnov test need not be an
# issue with large numbers of samples.Repeat the steps above, but for germination ratio as the demographic response variable
germ_trait_model_full <- glmmTMB(germ_ratio ~
ca_mg* log_sla +
depth*log_sla +
sand*log_sla +
ca_mg*log_srl +
depth*log_srl +
sand*log_srl +
ca_mg*log_maxht +
sand*log_maxht +
depth*log_maxht +
ca_mg*log_seedmass +
sand*log_seedmass +
depth*log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
weights = expected_germinants_l,
data = perf_env_trt,
family = binomial)
summary(germ_trait_model_full)## Family: binomial ( logit )
## Formula:
## germ_ratio ~ ca_mg * log_sla + depth * log_sla + sand * log_sla +
## ca_mg * log_srl + depth * log_srl + sand * log_srl + ca_mg *
## log_maxht + sand * log_maxht + depth * log_maxht + ca_mg *
## log_seedmass + sand * log_seedmass + depth * log_seedmass +
## (1 | sp_code) + (1 | plot_num)
## Zero inflation: ~1
## Data: perf_env_trt
## Weights: expected_germinants_l
##
## AIC BIC logLik deviance df.resid
## 15329.4 15449.3 -7641.7 15283.4 1337
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## sp_code (Intercept) 0.93011 0.9644
## plot_num (Intercept) 0.09286 0.3047
## Number of obs: 1360, groups: sp_code, 17; plot_num, 16
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.436585 0.246548 -5.827 5.65e-09 ***
## ca_mg 0.111468 0.092632 1.203 0.22885
## log_sla -0.291001 0.240677 -1.209 0.22663
## depth 0.019931 0.084204 0.237 0.81289
## sand -0.228518 0.085386 -2.676 0.00744 **
## log_srl -0.252324 0.256297 -0.984 0.32487
## log_maxht -0.316444 0.252178 -1.255 0.20954
## log_seedmass 0.578586 0.241472 2.396 0.01657 *
## ca_mg:log_sla 0.047782 0.010864 4.398 1.09e-05 ***
## log_sla:depth -0.013834 0.010175 -1.360 0.17395
## log_sla:sand 0.008868 0.011166 0.794 0.42707
## ca_mg:log_srl 0.082938 0.011100 7.472 7.91e-14 ***
## depth:log_srl -0.046089 0.010374 -4.443 8.87e-06 ***
## sand:log_srl -0.011107 0.011102 -1.000 0.31710
## ca_mg:log_maxht 0.086169 0.013022 6.617 3.67e-11 ***
## sand:log_maxht -0.029945 0.014394 -2.080 0.03749 *
## depth:log_maxht -0.062612 0.012562 -4.984 6.22e-07 ***
## ca_mg:log_seedmass 0.084725 0.012397 6.834 8.26e-12 ***
## sand:log_seedmass 0.029958 0.013294 2.253 0.02423 *
## depth:log_seedmass -0.016727 0.011571 -1.446 0.14829
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0817 0.1021 -20.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1germ_trait_model_null <- glmmTMB(germ_ratio ~
1 +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
weights = expected_germinants_l,
family = binomial)
germ_trait_model_additive <- glmmTMB(germ_ratio ~
ca_mg + sand + depth +
log_sla + log_srl + log_maxht + log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
weights = expected_germinants_l,
data = perf_env_trt,
family = binomial)
bbmle::AICctab(germ_trait_model_full, germ_trait_model_null, germ_trait_model_additive)## dAICc df
## germ_trait_model_full 0.0 23
## germ_trait_model_additive 397.8 11
## germ_trait_model_null 401.7 4resids <- DHARMa::simulateResiduals(germ_trait_model_full)
plot(resids)## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
# Looks some issue with uniformity - though the QQ looks fairly
# linear, and Florian Hartig's notes (e.g. https://github.com/florianhartig/DHARMa/issues/181)
# suggest that significance of Kolmogorov-Smirnov test need not be an
# issue with large numbers of samples.Analysis of Germ*Seed_production
germprod_trait_model_full <- glmmTMB(germ_prod ~
ca_mg* log_sla +
depth*log_sla +
sand*log_sla +
ca_mg*log_srl +
depth*log_srl +
sand*log_srl +
ca_mg*log_maxht +
sand*log_maxht +
depth*log_maxht +
ca_mg*log_seedmass +
sand*log_seedmass +
depth*log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
family = nbinom2(link = "log"))
summary(germprod_trait_model_full)## Family: nbinom2 ( log )
## Formula:
## germ_prod ~ ca_mg * log_sla + depth * log_sla + sand * log_sla +
## ca_mg * log_srl + depth * log_srl + sand * log_srl + ca_mg *
## log_maxht + sand * log_maxht + depth * log_maxht + ca_mg *
## log_seedmass + sand * log_seedmass + depth * log_seedmass +
## (1 | sp_code) + (1 | plot_num)
## Zero inflation: ~1
## Data: perf_env_trt
##
## AIC BIC logLik deviance df.resid
## 9788.4 9913.5 -4870.2 9740.4 1336
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## sp_code (Intercept) 3.5032 1.8717
## plot_num (Intercept) 0.5117 0.7153
## Number of obs: 1360, groups: sp_code, 17; plot_num, 16
##
## Dispersion parameter for nbinom2 family (): 0.699
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.951120 0.491284 6.007 1.89e-09 ***
## ca_mg 0.698442 0.223851 3.120 0.00181 **
## log_sla -0.607563 0.468589 -1.297 0.19478
## depth 0.351933 0.201849 1.744 0.08124 .
## sand -0.101619 0.204659 -0.497 0.61952
## log_srl -0.111318 0.499155 -0.223 0.82353
## log_maxht 0.158034 0.492304 0.321 0.74820
## log_seedmass 1.019855 0.470353 2.168 0.03014 *
## ca_mg:log_sla 0.256983 0.053085 4.841 1.29e-06 ***
## log_sla:depth -0.005919 0.042368 -0.140 0.88889
## log_sla:sand 0.054600 0.048987 1.115 0.26503
## ca_mg:log_srl 0.098004 0.066579 1.472 0.14102
## depth:log_srl -0.040143 0.053394 -0.752 0.45215
## sand:log_srl -0.118087 0.054105 -2.183 0.02907 *
## ca_mg:log_maxht 0.052355 0.062661 0.836 0.40342
## sand:log_maxht 0.089636 0.058791 1.525 0.12734
## depth:log_maxht 0.034176 0.057026 0.599 0.54897
## ca_mg:log_seedmass -0.110793 0.061638 -1.797 0.07226 .
## sand:log_seedmass -0.136593 0.061481 -2.222 0.02630 *
## depth:log_seedmass 0.033622 0.052208 0.644 0.51958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1965 0.1096 -10.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1germprod_trait_model_null <- glmmTMB(germ_prod ~
1 +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
family = nbinom2(link = "log"))
germprod_trait_model_additive <- glmmTMB(germ_prod ~
ca_mg + sand + depth +
log_sla + log_srl + log_maxht + log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt,
family = nbinom2(link = "log"))
bbmle::AICctab(germprod_trait_model_null, germprod_trait_model_full, germprod_trait_model_additive)## dAICc df
## germprod_trait_model_full 0.0 24
## germprod_trait_model_additive 29.3 12
## germprod_trait_model_null 36.1 5# Check VIF in the conditional component of the full model
performance::check_collinearity(germprod_trait_model_full, component = "conditional")# Check other model diagnostics using DHARMa
resids <- DHARMa::simulateResiduals(germprod_trait_model_full)
plot(resids)
Now, make plots
## plot model coefficients ##
model_terms_full <-
data.frame(confint(trait_model_full)) %>%
rownames_to_column() %>%
rename(term = rowname, lower_bound = X2.5..,
upper_bound = X97.5.., estimate = Estimate) %>%
filter(str_detect(string = term, pattern = "cond.")) %>%
filter(!str_detect(string = term, pattern = "Intercept")) %>%
mutate(term = str_replace(term, "cond.", ""),
order = c(19,16,18,17,15:6,4,5,3,1,2)) %>%
arrange(order) %>%
mutate(term_clean = rev(c("Soil Ca:Mg", "Soil depth", "Sand content",
"Specific leaf area (SLA)", "Specific root length (SRL)", "Max height", "Seed mass",
"Ca:Mg x SLA", "Depth x SLA", "Sand content x SLA",
"Ca:Mg x SRL", "Depth x SRL", "Sand content x SRL",
"Ca:Mg x Max height", "Depth x Max height", "Sand content x Max height",
"Ca:Mg x Seed mass", "Depth x Seed mass", "Sand content x Seed mass"))) %>%
arrange(-order) %>%
mutate(term_clean = fct_reorder(term_clean, -order),
fill_color = ifelse(upper_bound > 0 & lower_bound < 0, "NS", "S"))
germ_model_terms_full <-
data.frame(confint(germ_trait_model_full)) %>%
rownames_to_column() %>%
rename(term = rowname, lower_bound = X2.5..,
upper_bound = X97.5.., estimate = Estimate) %>%
filter(str_detect(string = term, pattern = "cond.")) %>%
filter(!str_detect(string = term, pattern = "Intercept")) %>%
mutate(term = str_replace(term, "cond.", ""),
order = c(19,16,18,17,15:6,4,5,3,1,2)) %>%
arrange(order) %>%
mutate(term_clean = rev(c("Soil Ca:Mg", "Soil depth", "Sand content",
"Specific leaf area (SLA)", "Specific root length (SRL)", "Max height", "Seed mass",
"Ca:Mg x SLA", "Depth x SLA", "Sand content x SLA",
"Ca:Mg x SRL", "Depth x SRL", "Sand content x SRL",
"Ca:Mg x Max height", "Depth x Max height", "Sand content x Max height",
"Ca:Mg x Seed mass", "Depth x Seed mass", "Sand content x Seed mass"))) %>%
arrange(-order) %>%
mutate(term_clean = fct_reorder(term_clean, -order),
fill_color = ifelse(upper_bound > 0 & lower_bound < 0, "NS", "S"))
germprod_model_terms_full <-
data.frame(confint(germprod_trait_model_full)) %>%
rownames_to_column() %>%
rename(term = rowname, lower_bound = X2.5..,
upper_bound = X97.5.., estimate = Estimate) %>%
filter(str_detect(string = term, pattern = "cond.")) %>%
filter(!str_detect(string = term, pattern = "Intercept")) %>%
mutate(term = str_replace(term, "cond.", ""),
order = c(19,16,18,17,15:6,4,5,3,1,2)) %>%
arrange(order) %>%
mutate(term_clean = rev(c("Soil Ca:Mg", "Soil depth", "Sand content",
"Specific leaf area (SLA)", "Specific root length (SRL)", "Max height", "Seed mass",
"Ca:Mg x SLA", "Depth x SLA", "Sand content x SLA",
"Ca:Mg x SRL", "Depth x SRL", "Sand content x SRL",
"Ca:Mg x Max height", "Depth x Max height", "Sand content x Max height",
"Ca:Mg x Seed mass", "Depth x Seed mass", "Sand content x Seed mass"))) %>%
arrange(order) %>%
mutate(term_clean = fct_reorder(term_clean, -order),
fill_color = ifelse(upper_bound > 0 & lower_bound < 0, "NS", "S"))Make Figure 2
## Compute the conditional and marginal R2
rsq_germprod <- performance::r2_nakagawa(germprod_trait_model_full)
rsq_prod <- performance::r2_nakagawa(trait_model_full)
rsq_germ <- performance::r2_nakagawa(germ_trait_model_full)
# Combine the model coefs dfs
combined_model_coefs <- left_join(model_terms_full,
germ_model_terms_full, by = c("term", "term_clean", "order"), )
all_combined_model_coefs <- left_join(combined_model_coefs, germprod_model_terms_full,
by = c("term", "term_clean", "order"), suffix = c("1","2"))
# Make the plot
(all_combined_model_coefs_plot <-
ggplot(all_combined_model_coefs) +
geom_errorbar(aes(x = order, ymin = lower_bound.x,
ymax = upper_bound.x), color = alpha("grey", 0)) +
geom_vline(xintercept = c(2,5,11,12), size = 8, color = alpha("#0072B2", 0.5)) +
geom_vline(xintercept = c(7), size = 8, color = alpha("#D55E00", 0.25)) +
# Germination ratio
geom_point(aes(x = order, y = estimate.y, fill = fill_color.y),
size = 2.35, shape = 21, stroke = 1.5) +
geom_errorbar(aes(x = order, ymin = lower_bound.y,
ymax = upper_bound.y), width = 0) +
geom_hline(yintercept = 1.5) +
# Seed production
geom_point(aes(x = order, y = estimate.x + 3, fill = fill_color.x),
size = 2.35, shape = 21, stroke = 1.5) +
geom_errorbar(aes(x = order, ymin = lower_bound.x + 3,
ymax = upper_bound.x + 3), width = 0) +
geom_hline(yintercept = c(0,3,6), linetype = "dashed") +
geom_hline(yintercept = 4.45) +
# Seed production x Germination ratio
geom_point(aes(x = order, y = estimate + 6, fill = fill_color),
size = 2.35, shape = 21, stroke = 1.5) +
geom_errorbar(aes(x = order, ymin = lower_bound + 6,
ymax = upper_bound + 6), width = 0) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous(breaks = c(-1, 0, 1, 2, 3, 4, 5, 6,7,8),
labels = c("-1",
paste0("0\n\n Marginal R2 = ",round(rsq_germ$R2_marginal,3),
"\nConditional R2 = ",round(rsq_germ$R2_conditional,3)), "1",
"-1",paste0("0\n\n Marginal R2 = ",round(rsq_prod$R2_marginal,3),
"\nConditional R2 = ",round(rsq_prod$R2_conditional,3)),"1",
"-1",paste0("0\n\n Marginal R2 = ",round(rsq_germprod$R2_marginal,3),
"\nConditional R2 = ",round(rsq_germprod$R2_conditional,3)),"1","2"
),
sec.axis = dup_axis(breaks = c(0, 3, 6),
labels = c("(A) Germination rate (g)\n",
"(B) Fecundity (F)\n",
"(C) g*F\n"),
name = "")) +
scale_x_continuous(breaks = combined_model_coefs$order, labels =combined_model_coefs$term_clean, expand = c(0,0)) +
scale_fill_manual(values = c("white", "grey60")) +
xlab("") +
ylab("") +
theme(axis.text = element_text(color = "black", size = 11),
axis.title = element_text(size = 12),
# axis.line.x.top = element_line(size = 0),
axis.text.x.top = element_text(face = "bold"),
axis.line.x.top = element_blank(),
legend.position = "NA") +
coord_flip()
)
# ggsave("figures/all-combined-model-coefficients.pdf", all_combined_model_coefs_plot, height = 7, width = 8)## Interaction plots ##
colpal <- c( "#E69F00", "#009E73")
get_pval <- function(model_out, mod_term) {
pval <- model_out %>% broom.mixed::tidy() %>%
filter(term == mod_term) %>%
dplyr::select(p.value) %>% unlist
pval <- ifelse(pval < .05,
paste0("\t p = ", scales::pvalue(pval)),
"\t N.S.")
pval
}
# make variables to help un-scale the axes
sla_sd <- sd(perf_env_trt$unscaled_log_sla)
sla_mean <- mean(perf_env_trt$unscaled_log_sla)
srl_sd <- sd(perf_env_trt$unscaled_log_srl)
srl_mean <- mean(perf_env_trt$unscaled_log_srl)
maxht_sd <- sd(perf_env_trt$unscaled_log_maxht)
maxht_mean <- mean(perf_env_trt$unscaled_log_maxht)
seedmass_sd <- sd(perf_env_trt$unscaled_log_seedmass)
seedmass_mean <- mean(perf_env_trt$unscaled_log_seedmass)
ca_mg_min <- round(min(perf_env_trt$unscaled_ca_mg), 3)
ca_mg_max <- round(max(perf_env_trt$unscaled_ca_mg), 3)
sand_min <- round(min(perf_env_trt$unscaled_sand), 3)
sand_max <- round(max(perf_env_trt$unscaled_sand), 3)
depth_min <- round(min(perf_env_trt$unscaled_depth), 3)
depth_max <- round(max(perf_env_trt$unscaled_depth), 3)
# FECUNDITY ONLY --------
# SLA:CA_MG
eff_sla_camg_df <- effects::effect("ca_mg:log_sla", trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_sla = c(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla)))) %>%
data.frame
# The warning message can be ignored according to Ben Bolker in the GLMMTMB docs:
# https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/inst/doc/model_evaluation.pdf
plot_sla_camg <-
ggplot(eff_sla_camg_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha=factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_sla", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_sla_depth_df <- effects::effect("log_sla:depth", trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_sla = c(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla)))) %>%
data.frame
plot_sla_depth <-
ggplot(eff_sla_depth_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower, ymax=upper,fill=factor(depth), alpha=factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "log_sla:depth", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_srl_sand_df <- effects::effect("sand:log_srl", trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_srl = c(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl)))) %>%
data.frame
plot_srl_sand <-
ggplot(eff_srl_sand_df, aes(x = (log_srl*srl_sd) + srl_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_srl", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific root length)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_maxht_depth_df <- effects::effect("depth:log_maxht", trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_maxht = c(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht)))) %>%
data.frame
plot_maxht_depth <-
ggplot(eff_maxht_depth_df, aes(x = (log_maxht*maxht_sd) + maxht_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_maxht", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Maximum height)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# Seed mass x Depth
eff_maxht_seedmass_df <- effects::effect("depth:log_seedmass", trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_maxht = c(min(perf_env_trt$log_seedmass),
max(perf_env_trt$log_seedmass)))) %>%
data.frame
plot_maxht_seedmass <-
ggplot(eff_maxht_seedmass_df, aes(x = (log_seedmass*seedmass_sd) + seedmass_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_seedmass", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Seed mass)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))Trait-environment interaction in germprod
# g*F ----------
germprod_eff_sla_camg_df <- effects::effect("ca_mg:log_sla", germprod_trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_sla = c(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla)))) %>%
data.frame## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are
## one-column matrices that were converted to vectorsgermprod_plot_sla_camg <-
ggplot(germprod_eff_sla_camg_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha=factor(ca_mg)), colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_sla", model_out = germprod_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("g*F") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
germprod_eff_sla_depth_df <- effects::effect("log_sla:depth", germprod_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_sla = c(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla)))) %>%
data.frame## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): the predictors
## ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are one-column
## matrices that were converted to vectorsgermprod_plot_sla_depth <-
ggplot(germprod_eff_sla_depth_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower, ymax=upper,fill=factor(depth), alpha=factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "log_sla:depth", model_out = germprod_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("g*F") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
germprod_eff_srl_sand_df <- effects::effect("sand:log_srl", germprod_trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_srl = c(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl)))) %>%
data.frame## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): the predictors
## ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are one-column
## matrices that were converted to vectorsgermprod_plot_srl_sand <-
ggplot(germprod_eff_srl_sand_df, aes(x = (log_srl*srl_sd) + srl_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(sand), alpha = factor(sand)), colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_srl", model_out = germprod_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific root length)") +
ylab("g*F") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
germprod_eff_maxht_depth_df <- effects::effect("depth:log_maxht", germprod_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_maxht = c(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht)))) %>%
data.frame## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): the predictors
## ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are one-column
## matrices that were converted to vectorsgermprod_plot_maxht_depth <-
ggplot(germprod_eff_maxht_depth_df, aes(x = (log_maxht*maxht_sd) + maxht_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(depth), alpha=factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_maxht", model_out = germprod_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Maximum height)") +
ylab("g*F") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# Seed mass x Depth
germprod_eff_depth_seedmass_df <- effects::effect("depth:log_seedmass", germprod_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_maxht = c(min(perf_env_trt$log_seedmass),
max(perf_env_trt$log_seedmass)))) %>%
data.frame## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): the predictors
## ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are one-column
## matrices that were converted to vectorsgermprod_plot_maxht_seedmass <-
ggplot(germprod_eff_depth_seedmass_df, aes(x = (log_seedmass*seedmass_sd) + seedmass_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(depth),alpha=factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_seedmass", model_out = germprod_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Seed mass) (g)") +
ylab("g*F") +
theme(legend.position = "top",
legend.text = element_text(size = 12))Trait-environment interaction in the germination ratios
# GERMINATION RATES ---------
# SLA v CA MG as predictor of seed germination
germ_eff_sla_camg_df <- effects::effect("ca_mg:log_sla", germ_trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_sla = seq(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are
## one-column matrices that were converted to vectorsgerm_plot_sla_camg <-
ggplot(germ_eff_sla_camg_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha=factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_sla", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# SLA vs Depth
germ_eff_sla_depth_df <- effects::effect("log_sla:depth", germ_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_sla = seq(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_sla_depth <-
ggplot(germ_eff_sla_depth_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "log_sla:depth", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# SRL v Sand
germ_eff_srl_sand_df <- effects::effect("sand:log_srl", germ_trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_srl = seq(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_srl_sand <-
ggplot(germ_eff_srl_sand_df, aes(x = (log_srl*srl_sd) + srl_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_srl", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Sand content = ", sand_min, "%"),
paste0("Sand content = ", sand_max, "%")),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Specific root length)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# Depth x Max Height
germ_eff_maxht_depth_df <- effects::effect("depth:log_maxht", germ_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_maxht = seq(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht),0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_maxht_depth <-
ggplot(germ_eff_maxht_depth_df, aes(x = (log_maxht*maxht_sd) + maxht_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_maxht", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Maximum height)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# Seed mass x Depth
germ_eff_maxht_seedmass_df <- effects::effect("depth:log_seedmass", germ_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_seedmass = seq(min(perf_env_trt$log_seedmass),
max(perf_env_trt$log_seedmass), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_maxht_seedmass <-
ggplot(germ_eff_maxht_seedmass_df, aes(x = (log_seedmass*seedmass_sd) + seedmass_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_seedmass", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min, "cm"),
paste0("Soil depth = ", depth_max, "cm")),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Seed mass)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))Combine all the plots above into one
lsla_camg_plot <- lsla_camg_plot + ylab(TeX("CWM log(SLA) (cm^2/mg)"))
lsla_depth_plot <- lsla_depth_plot + ylab(TeX("CWM log(SLA) (cm^2/mg)"))
lsrl_sand_plot <- lsrl_sand_plot + ylab("CWM log(SRL) (m/g)")
lmaxht_depth_plot <- lmaxht_depth_plot + ylab("CWM log(Max height) (cm)")
lseedmass_depth_plot <- lseedmass_depth_plot + ylab("CWM log(Seed mass) (g)")
fig3 <-
{lsla_camg_plot + lsla_depth_plot + lsrl_sand_plot + lmaxht_depth_plot + lseedmass_depth_plot +
germprod_plot_sla_camg + germprod_plot_sla_depth + germprod_plot_srl_sand +
germprod_plot_maxht_depth + germprod_plot_maxht_seedmass} +
plot_layout(ncol = 2, byrow = F) +
plot_annotation(tag_levels = "A")
# ggsave("figures/cwm-tXe-plots1.pdf", fig3, height = 12, width = 9)
supplemental_tXefigure <-
{(germ_plot_sla_camg + labs(tag = "A") + plot_sla_camg + labs(tag = "F")+ plot_layout(guides = "collect"))} /
{(germ_plot_sla_depth + labs(tag = "B") + plot_sla_depth + labs(tag = "G") + plot_layout(guides = "collect"))} /
{(germ_plot_srl_sand + labs(tag = "C") + plot_srl_sand + labs(tag = "H") + plot_layout(guides = "collect"))} /
{(germ_plot_maxht_depth + labs(tag = "D") + plot_maxht_depth + labs(tag = "I") + plot_layout(guides = "collect"))} /
{(germ_plot_maxht_seedmass + labs(tag = "E") + plot_maxht_seedmass + labs(tag = "J") + plot_layout(guides = "collect"))} +
plot_layout(byrow = T) &
theme(legend.position = "top",
plot.subtitle = element_text(face = "bold"))
# ggsave("figures/cwm-tXe-plots-germ-and-prod.pdf", supplemental_tXefigure, height = 12, width = 9)Now, make plots of all the non-significant T-E interactions:
# Make plots of all non-sig tXe ----
eff_sla_sand_df <- effects::effect("log_sla:sand", trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_sla = c(min(perf_env_trt$log_sla),
max(perf_env_trt$log_sla)))) %>%
data.frame
plot_sla_sand <-
ggplot(eff_sla_sand_df, aes(x = (log_sla*sla_sd)+sla_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "log_sla:sand", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil sand = ", sand_min),
paste0("Soil sand = ", sand_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil sand = ", sand_min),
paste0("Soil sand = ", sand_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil sand = ", sand_min),
paste0("Soil sand = ", sand_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific leaf area)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_srl_camg_df <- effects::effect("ca_mg:log_srl", trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_srl = c(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl)))) %>%
data.frame
plot_srl_camg <-
ggplot(eff_srl_camg_df, aes(x = (log_srl*srl_sd)+sla_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha = factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_srl", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific root length)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_srl_depth_df <- effects::effect("depth:log_srl", trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_srl = c(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl)))) %>%
data.frame
plot_srl_depth <-
ggplot(eff_srl_depth_df, aes(x = (log_srl*srl_sd)+sla_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_srl", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Depth = ", depth_min),
paste0("Soil Depth = ", depth_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Depth = ", depth_min),
paste0("Soil Depth = ", depth_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Depth = ", depth_min),
paste0("Soil Depth = ", depth_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Specific root length)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_maxht_sand_df <- effects::effect("sand:log_maxht", trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_maxht = c(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht)))) %>%
data.frame
plot_maxht_sand <-
ggplot(eff_maxht_sand_df, aes(x = (log_maxht*maxht_sd)+maxht_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_maxht", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil sand = ", sand_min),
paste0("Soil sand = ", sand_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil sand = ", sand_min),
paste0("Soil sand = ", sand_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil sand = ", sand_min),
paste0("Soil sand = ", sand_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Max height)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
eff_maxht_camg_df <- effects::effect("ca_mg:log_maxht", trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_maxht = c(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht)))) %>%
data.frame
plot_maxht_camg <-
ggplot(eff_maxht_camg_df, aes(x = (log_maxht*maxht_sd)+maxht_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha=factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_maxht", model_out = trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Max height)") +
ylab("Fecundity") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
comboplots_supplemental <-
{lsrl_camg_plot + lsrl_depth_plot + lsla_sand_plot +
lmaxht_sand_plot + lmaxht_camg_plot + plot_srl_camg +
plot_srl_depth + plot_sla_sand + plot_maxht_sand +
plot_maxht_camg} +
plot_layout(ncol = 2, byrow = F) +
plot_annotation(tag_levels = "A")
# ggsave("figures/supplemental/cwm-tXe-plots-supplemental.pdf", comboplots_supplemental, height = 12, width = 9)
comboplots_supplemental
Now, make plots of the TxE interactions that come out as significant in the vital rates analyses but not in the CWM traits:
# Interaction plots for significant interactions in the model
# (But no evidence at the CWM-level)
germprod_eff_seedmass_sand_df <- effects::effect("sand:log_seedmass", germprod_trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_maxht = c(min(perf_env_trt$log_seedmass),
max(perf_env_trt$log_seedmass)))) %>%
data.frame## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects: computed variances may be incorrect## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictors ca_mg, log_sla, depth, sand, log_srl, log_maxht, log_seedmass are
## one-column matrices that were converted to vectorsgermprod_plot_seedmass_sand <-
ggplot(germprod_eff_seedmass_sand_df, aes(x = (log_seedmass*seedmass_sd) + seedmass_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_seedmass", model_out = germprod_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil sand content = ", sand_min, "%"),
paste0("Soil sand content = ", sand_max, "%")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil sand content = ", sand_min, "%"),
paste0("Soil sand content = ", sand_max, "%")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil sand content = ", sand_min, "%"),
paste0("Soil sand content = ", sand_max, "%")),
values = c(0.25,0.5)) +
scale_y_log10() +
xlab("log(Seed mass) (g)") +
ylab("g*F") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# _---------
germ_eff_srl_depth_df <- effects::effect("depth:log_srl", germ_trait_model_full,
xlevels = list(depth = c(min(perf_env_trt$depth),
max(perf_env_trt$depth)),
log_srl = seq(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_srl_depth <-
ggplot(germ_eff_srl_depth_df, aes(x = (log_srl*srl_sd)+srl_mean, y = fit, color=factor(depth)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(depth), alpha = factor(depth)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "depth:log_srl", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min),
paste0("Soil depth = ", depth_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min),
paste0("Soil depth = ", depth_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil depth = ", depth_min),
paste0("Soil depth = ", depth_max)),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Specific root length)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# ------
germ_eff_srl_camg_df <- effects::effect("ca_mg:log_srl", germ_trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_srl = seq(min(perf_env_trt$log_srl),
max(perf_env_trt$log_srl), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_srl_camg <-
ggplot(germ_eff_srl_camg_df, aes(x = (log_srl*srl_sd)+srl_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha = factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_srl", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Specific root length)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# CaMg Max Ht -----------
germ_eff_maxht_camg_df <- effects::effect("ca_mg:log_maxht", germ_trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_maxht = seq(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_maxht_camg <-
ggplot(germ_eff_maxht_camg_df, aes(x = (log_maxht*maxht_sd)+maxht_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha = factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_maxht", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Max. Height)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# Sand Max Ht -----------
germ_eff_maxht_sand_df <- effects::effect("sand:log_maxht", germ_trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_maxht = seq(min(perf_env_trt$log_maxht),
max(perf_env_trt$log_maxht), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_maxht_sand <-
ggplot(germ_eff_maxht_sand_df, aes(x = (log_maxht*maxht_sd)+maxht_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_maxht", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Sand Content = ", sand_min, "%"),
paste0("Soil Sand Content = ", sand_max, "%")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Sand Content = ", sand_min, "%"),
paste0("Soil Sand Content = ", sand_max, "%")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Sand Content = ", sand_min, "%"),
paste0("Soil Sand Content = ", sand_max, "%")),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Max. Height)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# Sand content seed mass----
germ_eff_seedmass_sand_df <- effects::effect("sand:log_seedmass", germ_trait_model_full,
xlevels = list(sand = c(min(perf_env_trt$sand),
max(perf_env_trt$sand)),
log_seedmass = seq(min(perf_env_trt$log_seedmass),
max(perf_env_trt$log_seedmass), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_seedmass_sand <-
ggplot(germ_eff_seedmass_sand_df, aes(x = (log_seedmass*seedmass_sd)+seedmass_mean, y = fit, color=factor(sand)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(sand), alpha = factor(sand)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "sand:log_seedmass", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Sand Content = ", sand_min, "%"),
paste0("Soil Sand Content = ", sand_max, "%")),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Sand Content = ", sand_min, "%"),
paste0("Soil Sand Content = ", sand_max, "%")),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Sand Content = ", sand_min, "%"),
paste0("Soil Sand Content = ", sand_max, "%")),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Seed Mass)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# CaMg seed mass----
germ_eff_seedmass_camg_df <- effects::effect("ca_mg:log_seedmass", germ_trait_model_full,
xlevels = list(ca_mg = c(min(perf_env_trt$ca_mg),
max(perf_env_trt$ca_mg)),
log_seedmass = seq(min(perf_env_trt$log_seedmass),
max(perf_env_trt$log_seedmass), 0.1))) %>%
data.frame## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): the predictors ca_mg, log_sla, depth, sand,
## log_srl, log_maxht, log_seedmass are one-column matrices that were converted to
## vectorsgerm_plot_seedmass_camg <-
ggplot(germ_eff_seedmass_camg_df, aes(x = (log_seedmass*seedmass_sd)+seedmass_mean, y = fit, color=factor(ca_mg)))+
geom_line(size=1.25)+
geom_ribbon(aes(ymin=lower,ymax=upper, fill=factor(ca_mg), alpha = factor(ca_mg)),
colour= NA) +
annotate("text", x = -Inf, y = Inf, hjust = 0, vjust = 1,
label = get_pval(mod_term = "ca_mg:log_seedmass", model_out = germ_trait_model_full)) +
scale_fill_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_color_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = colpal) +
scale_alpha_manual(name = "",
labels = c(paste0("Soil Ca:Mg = ", ca_mg_min),
paste0("Soil Ca:Mg = ", ca_mg_max)),
values = c(0.25,0.5)) +
# scale_y_log10() +
xlab("log(Seed Mass)") +
ylab("Germination ratio") +
theme(legend.position = "top",
legend.text = element_text(size = 12))
# lsrl_camg_plot
lsrl_camg_plot <- lsrl_camg_plot + stat_smooth(method = "lm", formula = y ~ x + I(x^2))
# germprod_plot_seedmass_sand + lseedmass_sand_plot
# germ_plot_srl_depth + lsrl_depth_plot
# germ_plot_srl_camg + lsrl_camg_plot
# germ_plot_maxht_camg + lmaxht_camg_plot
# germ_plot_maxht_sand + lmaxht_sand_plot
# germ_plot_seedmass_sand + lseedmass_sand_plot
# germ_plot_seedmass_camg + lseedmass_camg_plot
comboplots_supplemental2 <-
{ germprod_plot_seedmass_sand + germ_plot_srl_depth + germ_plot_srl_camg +
germ_plot_maxht_camg + germ_plot_maxht_sand + germ_plot_seedmass_sand + germ_plot_seedmass_camg +
lseedmass_sand_plot + lsrl_depth_plot + lsrl_camg_plot +
lmaxht_camg_plot + lmaxht_sand_plot + lseedmass_sand_plot + lseedmass_camg_plot
} +
plot_layout(ncol = 2, byrow = F) +
plot_annotation(tag_levels = "A")
# ggsave("figures/supplemental/cwm-tXe-plots-supplemental2.pdf", comboplots_supplemental2, height = 15, width = 9)
comboplots_supplemental2
Supplemental analyses
Quadratic CWM trait-env models
p17_sitesum_env$ca_mg2 <- (p17_sitesum_env$ca_mg)^2
p17_sitesum_env$depth2 <- (p17_sitesum_env$depth)^2
p17_sitesum_env$sand2 <- (p17_sitesum_env$sand)^2
lsla_camg2_lm <- lm(lsla_cwm_mean~ca_mg + ca_mg2, data = p17_sitesum_env)
lsla_depth2_lm <- lm(lsla_cwm_mean~depth + depth2, data = p17_sitesum_env)
lsla_sand2_lm <- lm(lsla_cwm_mean~sand + sand2, data = p17_sitesum_env)
lsrl_camg2_lm <- lm(lsrl_cwm_mean~ca_mg + ca_mg2, data = p17_sitesum_env)
lsrl_depth2_lm <- lm(lsrl_cwm_mean~depth + depth2, data = p17_sitesum_env)
lsrl_sand2_lm <- lm(lsrl_cwm_mean~sand + sand2, data = p17_sitesum_env)
lmaxht_camg2_lm <- lm(lmaxht_cwm_mean~ca_mg + ca_mg2, data = p17_sitesum_env)
lmaxht_depth2_lm <- lm(lmaxht_cwm_mean~depth + depth2, data = p17_sitesum_env)
lmaxht_sand2_lm <- lm(lmaxht_cwm_mean~sand + sand2, data = p17_sitesum_env)
lseedmass_camg2_lm <- lm(lseedmass_cwm_mean~ca_mg + ca_mg2, data = p17_sitesum_env)
lseedmass_depth2_lm <- lm(lseedmass_cwm_mean~depth + depth2, data = p17_sitesum_env)
lseedmass_sand2_lm <- lm(lseedmass_cwm_mean~sand + sand2, data = p17_sitesum_env)
quadmodels_list <- list(lsla_camg2_lm, lsla_depth2_lm, lsla_sand2_lm,
lsrl_camg2_lm, lsrl_depth2_lm, lsrl_sand2_lm,
lmaxht_camg2_lm, lmaxht_depth2_lm, lmaxht_sand2_lm,
lseedmass_camg2_lm, lseedmass_depth2_lm, lseedmass_sand2_lm)
quadmodels_df <- map_df(quadmodels_list, tidy)
quadmodels_df$trait <- c("SLA", rep("", 8),
"SRL", rep("", 8),
"Max. Height", rep("", 8),
"Seed mass", rep("",8))
quadmodels_df$trait <- rep(c("SLA", "SRL", "Max. Height", "Seed mass"), each = 9)
quadmodels_df$environment <- rep(c("Ca:Mg" ,"","",
"Depth", "", "",
"Sand", "", ""), 4)
quadmodels_df$term <- rep(c("Intercept", "Linear", "Quadratic"), 12)
quadmodels_df$std.error <- round(quadmodels_df$std.error, 4)
quadmodels_df$statistic <- round(quadmodels_df$statistic, 4)
quadmodels_df$p.value <- round(quadmodels_df$p.value, 4)
quadmodels_df <- quadmodels_df %>%
dplyr::select(trait, environment, term, estimate, std.error, statistic,
p.value)
quadmodels_df_print <- knitr::kable(quadmodels_df, format = "latex", booktabs = T,
linesep = c(rep(c("", "", "\\addlinespace",
"", "", "\\addlinespace",
"", "", "\\addlinespace \\hline"), 3)),
caption = "Model output for quadratic relationships between CWM traits and environmental variables",
caption.short = "Testing for quadratic relationships between CWM traits and enironmental variables",
label = "CWMquadratic") %>%
kableExtra::kable_styling()
# saveRDS(object = quadmodels_df_print, file = "figures/quadmodels_df.Rds")Environment PCA
environment <- read_csv("data/all_environmental_data.csv") %>%
rename(plot_num = plot)## Rows: 16 Columns: 22## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): type, microsite
## dbl (20): plot, lat, lon, ele, Tmax, Tmin, organic_matter_ENR, pH, CEC_meq_1...##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.colnames(environment)## [1] "plot_num" "lat" "lon"
## [4] "ele" "type" "Tmax"
## [7] "Tmin" "organic_matter_ENR" "pH"
## [10] "CEC_meq_100g" "K_ppm" "Mg_ppm"
## [13] "Ca_ppm" "NH4_N_ppm" "Nitrate_ppm"
## [16] "soil_moisture" "sand" "clay"
## [19] "microsite" "depth" "P1_weak_bray_ppm"
## [22] "NaHCO3P_olsen_ppm"environment$ca_mg <- environment$Ca_ppm/environment$Mg_ppm
environment_pcacols <- environment %>%
dplyr::select(plot_num, organic_matter_ENR, pH, CEC_meq_100g,
K_ppm, Ca_ppm, Mg_ppm, NH4_N_ppm, Nitrate_ppm,
soil_moisture, sand, depth, Tmax, ca_mg, P1_weak_bray_ppm) %>%
mutate(plot_num = as.character(plot_num)) %>%
# scale the columns
mutate_if(is.numeric, scale) %>%
# plot number should be the name of the rows.
column_to_rownames("plot_num")
env_pca <- PCA(environment_pcacols, graph = F)
# Verify number of coordinates to retain using Parallel Analysis
paran::paran(environment_pcacols)##
## Using eigendecomposition of correlation matrix.
## Computing: 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
##
##
## Results of Horn's Parallel Analysis for component retention
## 420 iterations, using the mean estimate
##
## --------------------------------------------------
## Component Adjusted Unadjusted Estimated
## Eigenvalue Eigenvalue Bias
## --------------------------------------------------
## 1 3.170116 5.207093 2.036976
## 2 1.027142 2.447081 1.419939
## --------------------------------------------------
##
## Adjusted eigenvalues > 1 indicate dimensions to retain.
## (2 components retained)env_pca_coords <- env_pca$ind$coord
env_pca_coords <- env_pca_coords %>% data.frame %>% rownames_to_column("plot_num") %>%
mutate(plot_num = as.numeric(plot_num)) %>%
left_join(., environment %>%
dplyr::select(plot_num, microsite))## Joining, by = "plot_num"env_pca_12 <- factoextra::fviz_pca_biplot(env_pca, col.ind = "transparent") +
geom_point(data = env_pca_coords, aes(x = Dim.1, y = Dim.2, shape = microsite), size = 4)+
scale_shape_manual(values = c(20,21), labels = c("hummock", "matrix")) +
theme(legend.position = "bottom")+
labs(tile = "")
env_pca_12
# ggsave("figures/supplemental/env-pca.pdf", env_pca_12, height = 6, width = 6)Trait PCA
na_count <- sapply(strts, function(y) sum(length(which(is.na(y)))))
na_count## USDA_symbol leaf_size SLA LDMC seed_mass
## 0 1 1 1 1
## max_height SRL relative_spread phenology foliar_N
## 0 16 12 22 7
## CN_ratio d13C d15N dataset
## 7 7 7 0# Exclude traits for which we are lacking data from most
# species (leaf pH, TLP, rooting depth, seed 'size)
# also exclude notes/non-numeric columns for PCA
strts <- strts %>% dplyr::select(-dataset)
# dplyr::select(-c(leaf_pH, turgor_loss_point, dataset,
# max_height_data_source, LAI, LAR,
# seed_mass_data_source, notes, rooting_depth,
# seed_size, projected_area))
# Species codes for the focal species of the experiment
focal_sp_codes <- c("URLI5", "SACO6", "PLER3", "MICA",
"MEPO3", "LOWR2", "LACA7", "HOMU",
"HECO7","FEMI2", 'EUSP', "CLPU2",
'CLBO', 'CHGL', 'CEME2', 'BRMA3', 'AMME')
dominant_sp_codes <- c("AVFA", "BRDI3", "AVBA", "MIDO", "LOPE")
focal_sp_rownums <- which(strts$USDA_symbol %in% focal_sp_codes)
# Scale the numeric trait columns
strts_scaled <- strts %>%
column_to_rownames("USDA_symbol") %>%
replace_na(.,as.list(colMeans(.,na.rm=T))) %>%
mutate_if(is.numeric, scale)
trait_pca <- PCA(strts_scaled, graph = F)
focal_sp_coords <- data.frame(trait_pca$ind$coord[focal_sp_rownums,])
focal_sp_codes_sedgnames <- c("URLI", "SACO", "PLER", "MICA",
"MEPO", "ACWR", 'LACA', 'HOMU',
'HECO', 'VUMI', 'EUSP', 'CLPU',
'CLBO', 'CHGL', 'CEME', 'BRMA', 'AMME')
focal_sp_coords$sp_code <- focal_sp_codes_sedgnames
trait_pca_12 <- fviz_pca_biplot(trait_pca, col.ind = "grey", label = "var",
col.var = c("steelblue", "red", "steelblue", "red", "red", "red",
rep("steelblue", 6))) +
labs(title = "") +
geom_point(data = focal_sp_coords, aes(x = Dim.1, y = Dim.2),
color = "darkred", size = 2.5) +
ggrepel::geom_text_repel(data = focal_sp_coords, aes(x = Dim.1, y = Dim.2, label = sp_code)) +
theme(legend.position = "none")
trait_pca_12
# ggsave("figures/supplemental/trait-pca.pdf",trait_pca_12, height = 6, width = 6)
maxht_hist <- ggplot(strts) +
geom_histogram(aes(x = max_height), bins = 20) +
geom_rug(aes(x = max_height), color = "#E69F00", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% focal_sp_codes)) +
geom_rug(aes(x = max_height), color = "#0072B2", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% dominant_sp_codes)) +
scale_x_log10() +
xlab("Maximum height (cm)")
sla_hist <- ggplot(strts) +
geom_histogram(aes(x = SLA), bins = 20) +
geom_rug(aes(x = SLA), color = "#E69F00", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% focal_sp_codes)) +
geom_rug(aes(x = SLA), color = "#0072B2", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% dominant_sp_codes)) +
scale_x_log10() +
xlab(TeX("Specific leaf area (cm^2/g)"))
srl_hist <- ggplot(strts) +
geom_histogram(aes(x = SRL), bins = 20) +
geom_rug(aes(x = SRL), color = "#E69F00", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% focal_sp_codes)) +
geom_rug(aes(x = SRL), color = "#0072B2", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% dominant_sp_codes)) +
scale_x_log10() +
xlab("Specific root length (m/g)")
seedmass_hist <- ggplot(strts) +
geom_histogram(aes(x = seed_mass), bins = 20) +
geom_rug(aes(x = seed_mass), color = "#E69F00", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% focal_sp_codes)) +
geom_rug(aes(x = seed_mass), color = "#0072B2", size = 2, alpha = .75,
data = strts %>% filter(USDA_symbol %in% dominant_sp_codes)) +
scale_x_log10() +
xlab("Seed mass (g)")
fig_traithistograms <- {maxht_hist + sla_hist + srl_hist + seedmass_hist} +
plot_annotation(tag_levels = "A")
fig_traithistograms## Warning: Removed 1 rows containing non-finite values (stat_bin).## Warning: Removed 16 rows containing non-finite values (stat_bin).## Warning: Removed 1 rows containing non-finite values (stat_bin).
# ggsave("figures/supplemental/trait-histograms.pdf", fig_traithistograms, width = 7, height = 7)
# Calculate trait ranges
focal_sp_trts <- strts %>% filter(USDA_symbol %in% focal_sp_codes)
diff(range(log(focal_sp_trts$SRL), na.rm = T))/diff(range(log(strts$SRL), na.rm = T))## [1] 0.4373217diff(range(log(focal_sp_trts$SLA), na.rm = T))/diff(range(log(strts$SLA), na.rm = T))## [1] 0.6517279diff(range(log(focal_sp_trts$max_height), na.rm = T))/diff(range(log(strts$max_height), na.rm = T))## [1] 0.6324735diff(range(log(focal_sp_trts$seed_mass), na.rm = T))/diff(range(log(strts$seed_mass), na.rm = T))## [1] 0.7517895trait_table <- data.frame(
Organ = c("Whole plant"," "," "," ", " ", "Leaf", "", "", "", "", "Root"),
Trait = c("Max. height", "Canopy shape index", "Carbon isotope composition (dC13)",
"Phenology", "Seed mass", "Leaf size", "Specific leaf area",
"Leaf dry matter content", "C:N ratio", "Leaf N concentration",
"Specific root length"),
Units = c("cm", "dimensionless", "dC13", "day of year", "mg", "cm$^2$",
"g/cm$^2$", "mg/g", "dimensionless", "mg/g", "m/g")
)
trait_table_print <- knitr::kable(trait_table, format = "latex", booktabs = T,
label = "traittable",
caption = "List of traits measured for this study") %>%
kableExtra::kable_styling()
# saveRDS(object = trait_table_print, file = "figures/trait_table.Rds")Trait relationships with environmental PC axes
# CWM relationship with PCs 1 and 2 --------------
env_pca_coords$plot_num <- as.factor(env_pca_coords$plot_num)
p17_sitesum_env$plot_num <- as.factor(as.character(p17_sitesum_env$site))
p17_sitesum_env <- left_join(p17_sitesum_env, env_pca_coords, by = "plot_num")
lsla_dim1_lm <- make_trait_env_lms(trait = "lsla_cwm_mean", env = "Dim.1")
lsla_dim1_plot <- make_trait_env_plot(trait = "lsla_cwm_mean", env = "Dim.1",
lsla_dim1_lm, xl = "Env PCA Dim1", yl = "CWM log(SLA (cm^2/g))")
lsla_dim2_lm <- make_trait_env_lms(trait = "lsla_cwm_mean", env = "Dim.2")
lsla_dim2_plot <- make_trait_env_plot(trait = "lsla_cwm_mean", env = "Dim.2",
lsla_dim2_lm, xl = "Env PCA Dim2", yl = "CWM log(SLA (cm^2/g))")
lsrl_dim1_lm <- make_trait_env_lms(trait = "lsrl_cwm_mean", env = "Dim.1")
lsrl_dim1_plot <- make_trait_env_plot(trait = "lsrl_cwm_mean", env = "Dim.1",
lsrl_dim1_lm, xl = "Env PCA Dim1", yl = "CWM log(SRL (m/g))")
lsrl_dim2_lm <- make_trait_env_lms(trait = "lsrl_cwm_mean", env = "Dim.2")
lsrl_dim2_plot <- make_trait_env_plot(trait = "lsrl_cwm_mean", env = "Dim.2",
lsrl_dim2_lm, xl = "Env PCA Dim2", yl = "CWM log(SRL (m/g))")
lmaxht_dim1_lm <- make_trait_env_lms(trait = "lmaxht_cwm_mean", env = "Dim.1")
lmaxht_dim1_plot <- make_trait_env_plot(trait = "lmaxht_cwm_mean", env = "Dim.1",
lmaxht_dim1_lm, xl = "Env PCA Dim1", yl = "CWM log(Max Ht (cm))")
lmaxht_dim2_lm <- make_trait_env_lms(trait = "lmaxht_cwm_mean", env = "Dim.2")
lmaxht_dim2_plot <- make_trait_env_plot(trait = "lmaxht_cwm_mean", env = "Dim.2",
lmaxht_dim2_lm, xl = "Env PCA Dim2", yl = "CWM log(Max Ht (cm))")
lseedmass_dim1_lm <- make_trait_env_lms(trait = "lseedmass_cwm_mean", env = "Dim.1")
lseedmass_dim1_plot <- make_trait_env_plot(trait = "lseedmass_cwm_mean", env = "Dim.1",
lseedmass_dim1_lm, xl = "Env PCA Dim1", yl = "CWM log(Seed mass (mg))")
lseedmass_dim2_lm <- make_trait_env_lms(trait = "lseedmass_cwm_mean", env = "Dim.2")
lseedmass_dim2_plot <- make_trait_env_plot(trait = "lseedmass_cwm_mean", env = "Dim.2",
lseedmass_dim2_lm, xl = "Env PCA Dim2", yl = "CWM log(Seed mass (mg))")
trait_variation_PCs <-
lsla_dim1_plot + lsla_dim2_plot +
lsrl_dim1_plot + lsrl_dim2_plot +
lmaxht_dim1_plot + lmaxht_dim2_plot +
lseedmass_dim1_plot + lseedmass_dim2_plot +
plot_layout(ncol = 2)
trait_variation_PCs
# ggsave("figures/supplemental/trait-variation-PCs.pdf", trait_variation_PCs,height = 8, width = 6)
# Trait X PC interactions in GLMM
perf_env_trt_pca <- left_join(perf_env_trt, env_pca_coords, by = "plot_num")
trait_model_full_pca <- glmmTMB(num_seeds_produced ~
Dim.1* log_sla +
Dim.2*log_sla +
Dim.1*log_srl +
Dim.2*log_srl +
Dim.1*log_maxht +
Dim.2*log_maxht +
Dim.1*log_seedmass +
Dim.2*log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt_pca,
family = nbinom2(link = "log"))
summary(trait_model_full_pca)## Family: nbinom2 ( log )
## Formula:
## num_seeds_produced ~ Dim.1 * log_sla + Dim.2 * log_sla + Dim.1 *
## log_srl + Dim.2 * log_srl + Dim.1 * log_maxht + Dim.2 * log_maxht +
## Dim.1 * log_seedmass + Dim.2 * log_seedmass + (1 | sp_code) +
## (1 | plot_num)
## Zero inflation: ~1
## Data: perf_env_trt_pca
##
## AIC BIC logLik deviance df.resid
## 13248.1 13347.7 -6605.0 13210.1 1380
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## sp_code (Intercept) 0.4691 0.6849
## plot_num (Intercept) 0.1850 0.4301
## Number of obs: 1399, groups: sp_code, 17; plot_num, 16
##
## Dispersion parameter for nbinom2 family (): 1.51
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.25082 0.20236 25.948 < 2e-16 ***
## Dim.1 0.23705 0.04923 4.815 1.47e-06 ***
## log_sla -0.20603 0.17366 -1.186 0.235463
## Dim.2 0.20423 0.07148 2.857 0.004277 **
## log_srl 0.06014 0.18621 0.323 0.746698
## log_maxht 0.62917 0.18570 3.388 0.000704 ***
## log_seedmass 0.24108 0.17721 1.360 0.173692
## Dim.1:log_sla 0.03996 0.01452 2.753 0.005908 **
## log_sla:Dim.2 -0.01974 0.01833 -1.077 0.281671
## Dim.1:log_srl 0.02689 0.01613 1.667 0.095425 .
## Dim.2:log_srl 0.03199 0.02150 1.488 0.136814
## Dim.1:log_maxht -0.02524 0.01642 -1.537 0.124310
## Dim.2:log_maxht 0.04787 0.02376 2.015 0.043917 *
## Dim.1:log_seedmass 0.05053 0.01551 3.258 0.001121 **
## Dim.2:log_seedmass 0.03076 0.02179 1.411 0.158129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5343 0.0558 -9.577 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1germ_trait_model_full_pca <- glmmTMB(germ_ratio ~
Dim.1* log_sla +
Dim.2*log_sla +
Dim.1*log_srl +
Dim.2*log_srl +
Dim.1*log_maxht +
Dim.2*log_maxht +
Dim.1*log_seedmass +
Dim.2*log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
weights = expected_germinants_l,
data = perf_env_trt_pca,
family = binomial)
summary(germ_trait_model_full_pca)## Family: binomial ( logit )
## Formula:
## germ_ratio ~ Dim.1 * log_sla + Dim.2 * log_sla + Dim.1 * log_srl +
## Dim.2 * log_srl + Dim.1 * log_maxht + Dim.2 * log_maxht +
## Dim.1 * log_seedmass + Dim.2 * log_seedmass + (1 | sp_code) +
## (1 | plot_num)
## Zero inflation: ~1
## Data: perf_env_trt_pca
## Weights: expected_germinants_l
##
## AIC BIC logLik deviance df.resid
## 15319.2 15413.0 -7641.6 15283.2 1342
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## sp_code (Intercept) 0.9326 0.9657
## plot_num (Intercept) 0.1003 0.3168
## Number of obs: 1360, groups: sp_code, 17; plot_num, 16
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.459752 0.247793 -5.891 3.84e-09 ***
## Dim.1 0.129183 0.035055 3.685 0.000229 ***
## log_sla -0.294034 0.240994 -1.220 0.222432
## Dim.2 -0.143381 0.051106 -2.806 0.005023 **
## log_srl -0.254103 0.256629 -0.990 0.322097
## log_maxht -0.335655 0.252528 -1.329 0.183788
## log_seedmass 0.578184 0.241785 2.391 0.016788 *
## Dim.1:log_sla 0.014678 0.004304 3.410 0.000650 ***
## log_sla:Dim.2 -0.007533 0.006202 -1.215 0.224498
## Dim.1:log_srl 0.034664 0.004278 8.102 5.41e-16 ***
## Dim.2:log_srl -0.032410 0.006185 -5.240 1.60e-07 ***
## Dim.1:log_maxht 0.044544 0.005605 7.947 1.92e-15 ***
## Dim.2:log_maxht -0.049387 0.007853 -6.289 3.20e-10 ***
## Dim.1:log_seedmass 0.026236 0.004942 5.309 1.10e-07 ***
## Dim.2:log_seedmass -0.005971 0.006987 -0.855 0.392767
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1488 0.1062 -20.24 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1germprod_trait_model_full_pca <- glmmTMB(germ_prod ~
Dim.1* log_sla +
Dim.2*log_sla +
Dim.1*log_srl +
Dim.2*log_srl +
Dim.1*log_maxht +
Dim.2*log_maxht +
Dim.1*log_seedmass +
Dim.2*log_seedmass +
(1|sp_code) + (1|plot_num),
ziformula = ~1,
data = perf_env_trt_pca,
family = nbinom2(link = "log"))
summary(germprod_trait_model_full_pca)## Family: nbinom2 ( log )
## Formula:
## germ_prod ~ Dim.1 * log_sla + Dim.2 * log_sla + Dim.1 * log_srl +
## Dim.2 * log_srl + Dim.1 * log_maxht + Dim.2 * log_maxht +
## Dim.1 * log_seedmass + Dim.2 * log_seedmass + (1 | sp_code) +
## (1 | plot_num)
## Zero inflation: ~1
## Data: perf_env_trt_pca
##
## AIC BIC logLik deviance df.resid
## 9789.7 9888.8 -4875.9 9751.7 1341
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## sp_code (Intercept) 3.3746 1.8370
## plot_num (Intercept) 0.5394 0.7344
## Number of obs: 1360, groups: sp_code, 17; plot_num, 16
##
## Dispersion parameter for nbinom2 family (): 0.714
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.993617 0.485160 6.170 6.81e-10 ***
## Dim.1 0.395984 0.083313 4.753 2.00e-06 ***
## log_sla -0.599197 0.459945 -1.303 0.192657
## Dim.2 0.149949 0.120479 1.245 0.213276
## log_srl -0.143236 0.489852 -0.292 0.769975
## log_maxht 0.222412 0.483342 0.460 0.645406
## log_seedmass 0.985061 0.461572 2.134 0.032831 *
## Dim.1:log_sla 0.090777 0.020767 4.371 1.24e-05 ***
## log_sla:Dim.2 -0.015231 0.025755 -0.591 0.554256
## Dim.1:log_srl 0.077773 0.022600 3.441 0.000579 ***
## Dim.2:log_srl 0.021943 0.030478 0.720 0.471551
## Dim.1:log_maxht -0.018647 0.023523 -0.793 0.427954
## Dim.2:log_maxht 0.054810 0.032078 1.709 0.087511 .
## Dim.1:log_seedmass -0.006538 0.022175 -0.295 0.768109
## Dim.2:log_seedmass -0.013611 0.031109 -0.438 0.661737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1473 0.1048 -10.95 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Make plots for the PC x trait interaction models
prod_model_terms_full_pca <-
data.frame(confint(trait_model_full_pca)) %>%
rownames_to_column() %>%
rename(term = rowname, lower_bound = X2.5..,
upper_bound = X97.5.., estimate = Estimate) %>%
filter(str_detect(string = term, pattern = "cond.")) %>%
filter(!str_detect(string = term, pattern = "Intercept")) %>%
mutate(term = str_replace(term, "cond.", ""),
order = c(14,12,13,11,10,9:1)) %>%
arrange(order) %>%
mutate(term_clean = c("Dim2 x Seed mass", "Dim1 x Seed mass",
"Dim2 x Max height", "Dim1 x Max height",
"Dim2 x SRL", "Dim1 x SRL",
"Dim2 x SLA", "Dim1 x SLA",
"Seed mass", "Max height", "SRL", "SLA",
"Dim2", "Dim1")) %>%
mutate(term_clean = fct_reorder(term_clean, -order),
fill_color = ifelse(upper_bound > 0 & lower_bound < 0, "NS", "S"))
germ_model_terms_full_pca <-
data.frame(confint(germ_trait_model_full_pca)) %>%
rownames_to_column() %>%
rename(term = rowname, lower_bound = X2.5..,
upper_bound = X97.5.., estimate = Estimate) %>%
filter(str_detect(string = term, pattern = "cond.")) %>%
filter(!str_detect(string = term, pattern = "Intercept")) %>%
mutate(term = str_replace(term, "cond.", ""),
order = c(14,12,13,11,10,9:1)) %>%
arrange(order) %>%
mutate(term_clean = c("Dim2 x Seed mass", "Dim1 x Seed mass",
"Dim2 x Max height", "Dim1 x Max height",
"Dim2 x SRL", "Dim1 x SRL",
"Dim2 x SLA", "Dim1 x SLA",
"Seed mass", "Max height", "SRL", "SLA",
"Dim2", "Dim1")) %>%
mutate(term_clean = fct_reorder(term_clean, -order),
fill_color = ifelse(upper_bound > 0 & lower_bound < 0, "NS", "S"))
germprod_model_terms_full_pca <-
data.frame(confint(germprod_trait_model_full_pca)) %>%
rownames_to_column() %>%
rename(term = rowname, lower_bound = X2.5..,
upper_bound = X97.5.., estimate = Estimate) %>%
filter(str_detect(string = term, pattern = "cond.")) %>%
filter(!str_detect(string = term, pattern = "Intercept")) %>%
mutate(term = str_replace(term, "cond.", ""),
order = c(14,12,13,11,10,9:1)) %>%
arrange(order) %>%
mutate(term_clean = c("Dim2 x Seed mass", "Dim1 x Seed mass",
"Dim2 x Max height", "Dim1 x Max height",
"Dim2 x SRL", "Dim1 x SRL",
"Dim2 x SLA", "Dim1 x SLA",
"Seed mass", "Max height", "SRL", "SLA",
"Dim2", "Dim1")) %>%
mutate(term_clean = fct_reorder(term_clean, -order),
fill_color = ifelse(upper_bound > 0 & lower_bound < 0, "NS", "S"))
# Combine the model coefs dfs
pca_combined_model_coefs <- left_join(prod_model_terms_full_pca,
germ_model_terms_full_pca, by = c("term", "term_clean", "order"), )
pca_all_combined_model_coefs <- left_join(pca_combined_model_coefs, germprod_model_terms_full_pca,
by = c("term", "term_clean", "order"), suffix = c("1","2"))
rsq_germprod_pca <- performance::r2_nakagawa(germprod_trait_model_full_pca)
rsq_prod_pca <- performance::r2_nakagawa(trait_model_full_pca)
rsq_germ_pca <- performance::r2_nakagawa(germ_trait_model_full_pca)
# Make the plot
(pca_all_combined_model_coefs_plot <-
ggplot(pca_all_combined_model_coefs) +
geom_errorbar(aes(x = order, ymin = lower_bound.x,
ymax = upper_bound.x), color = alpha("grey", 0)) +
geom_vline(xintercept = c(3,4,5,8), size = 8, color = alpha("#0072B2", 0.5)) +
geom_vline(xintercept = c(2), size = 8, color = alpha("#D55E00", 0.25)) +
# Germination ratio
geom_point(aes(x = order, y = estimate.y, fill = fill_color.y),
size = 2.35, shape = 21, stroke = 1.5) +
geom_errorbar(aes(x = order, ymin = lower_bound.y,
ymax = upper_bound.y), width = 0) +
geom_hline(yintercept = 1.5) +
# Seed production
geom_point(aes(x = order, y = estimate.x + 3, fill = fill_color.x),
size = 2.35, shape = 21, stroke = 1.5) +
geom_errorbar(aes(x = order, ymin = lower_bound.x + 3,
ymax = upper_bound.x + 3), width = 0) +
geom_hline(yintercept = c(0,3,6), linetype = "dashed") +
geom_hline(yintercept = 4.45) +
# Seed production x Germination ratio
geom_point(aes(x = order, y = estimate + 6, fill = fill_color),
size = 2.35, shape = 21, stroke = 1.5) +
geom_errorbar(aes(x = order, ymin = lower_bound + 6,
ymax = upper_bound + 6), width = 0) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous(breaks = c(-1, 0, 1, 2, 3, 4, 5, 6,7,8),
# labels = c(-1,0,1, -1,0,1, -1,0,1,2),
labels = c("-1",
paste0("0\n\n Marginal R2 = ",round(rsq_germ_pca$R2_marginal,3),
"\nConditional R2 = ",round(rsq_germ_pca$R2_conditional,3)), "1",
"-1",paste0("0\n\n Marginal R2 = ",round(rsq_prod_pca$R2_marginal,3),
"\nConditional R2 = ",round(rsq_prod_pca$R2_conditional,3)),"1",
"-1",paste0("0\n\n Marginal R2 = ",round(rsq_germprod_pca$R2_marginal,3),
"\nConditional R2 = ",round(rsq_germprod_pca$R2_conditional,3)),"1","2"
),
sec.axis = dup_axis(breaks = c(0, 3, 6),
labels = c("(A) Germination rate (g)\n",
"(B) Fecundity (F)\n",
"(C) g*F\n"),
name = "")) +
scale_x_continuous(breaks = pca_all_combined_model_coefs$order, labels =pca_all_combined_model_coefs$term_clean, expand = c(0,0)) +
scale_fill_manual(values = c("white", "grey60")) +
xlab("") +
ylab("") +
theme(axis.text = element_text(color = "black", size = 11),
axis.title = element_text(size = 12),
# axis.line.x.top = element_line(size = 0),
axis.text.x.top = element_text(face = "bold"),
axis.line.x.top = element_blank(),
legend.position = "NA") +
coord_flip()
)
# ggsave("figures/supplemental/effect-sizes-pca.pdf", pca_all_combined_model_coefs_plot, height = 5, width = 8)